FINAL REPORT 


FOR 

NAG2-175 


\ — I 

ATMOSPHERE AND WIND MODELING FOR ATC 


PRINCIPAL INVESTIGATOR: 

Gary L. Slater 
Dept . AsE/EM 

University of Cincinnati 
Cincinnati, OH 45221 


November 1990 


(NASA-CR-l 96786) 
WINO MODELING FOR 
(C i nc i nnat i Un i v. ) 


ATMOSPHERIC AND 
ATC Final Report 
48 p 


N95-13725 
Unc 1 as 


63/ 47 0022282 





0 \ \& 



Atmospheric and Wind Modeling for ATC 


G.L, Slater 

University of Cincinnati 
Cincinnati, Ohio 


November 1990 


1. Atmospheric Modeling 


1.1 The standard atmosphere 

1.2 Atmosphere variations 

1.3 Atmosphere requirements for ATC 

1.4 Implementation of a software model for CTAS 

2. Wind modeling 

2.1 Wind data: NOAA Profiler system 

2.2 Wind profile estimation 

2.3 Incorporation of various data types into filtering scheme 

2.4 Spatial & temporal variation 

2.5 Software implementation into CTAS 

Appendix A: Mat lab code for atmospheric routines 
Appendix B: Matlab code for wind estimation 


1. Atmospheric Modeling 


Pressure, temperature and density of the atmosphere vary considerably 
with altitude, and in addition, exhibit additional spatial and temporal 
variations. Since aircraft motions and sensors are highly dependent on 
these atmospheric properties, it is important to understand and to model 
these variations in the development of an ATC advisory scheme. In this 
section we review briefly the concept of the model atmosphere, and how this 
model should be corrected to allow for real variations. 

1.1 The Model atmosphere 

To predict the gross pressure and density properties of the atmosphere 
at typical aircraft altitudes, It is sufficient to model the atmosphere as a 
perfect gas, and as a fluid in equilibrium on a non-rotating, flat earth. 
Such a fluid satisfies: 

(1) Perfect gas law: p = pRT 

(2) The hydrostatic equation: dp - -pgdh 

where p, p, T are the pressure, density and temperature respectively, R is 

the gas constant (Universal gas constant R divided by the molecular weight 
of air) and h is the geometric altitude measured positive upwards. 
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To determine the atmosphere completely, one additional relation between 
these variables is required. In a theoretical sense this should be a 
thermodynamic equation expressing conservation of energy; that is the heat 
balance between internal energy, and heat flux through radiation, conduction 
and convection. This last relationship is difficult to express; Instead, 
for the purposes of defining an "average” atmosphere, an average temperature 
variation with altitude is assumed, and this together with equations (l)-(2) 
define the model atmosphere. 

Various atmospheric models have been proposed and are in use. The most 
common and universally used is the ICAO (International Civil Aviation 
Authority) which was introduced in 1952 to cover altitudes up to 20km 
(approximately 65,500ft). The ICAO atmosphere Is similar to, but 
supercedes, the NACA (National Advisory Committee for Aeronautics) 
atmosphere developed In the United States and the ICAN (International 
Commission for Aerial Navigation) atmosphere developed In Europe. Since 
1952 this atmosphere has been refined to cover altitudes up to satellite 
levels but the levels up to 20km are essentially unchanged. (See, for 
example, the US Standard Atmosphere , 1962 [1]) This atmosphere is a 
statistical average over the year for mid- latitudes in the northern 
hemisphere . 

In the altitude range of interest the temperature is defined by two 
linear functions which define two atmospheric regions: 

a) Troposhere: 
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0 < h < 11km T - 288.15 - 0.0065h (deg K) 


b) Stratosphere: 

11km < h T— 216.65 (deg K) 

The mean temperature at the earth surface is 288.15 K - 15 C - 59 F - 
518.67 R . The slope of the temperature versus altitude curve is referred to 
as the lapse rate and in the troposphere has the defined value -0.0065 deg 
C/m ( or -0.00367 deg F/ft) . The stratosphere is modeled as an isothermal 
region at a temperature of -56.5 deg C (216.65 K) or -69.7 deg F (389.97 deg 
R) . The dividing line between the troposphere and the stratosphere is 
referred to as the tropopause and in the model atmosphere is defined at 11km 
— 36089ft (commonly rounded off to 36,000ft). The pressure at sea level 
(h-=0) is defined in the model atmosphere to be exactly p 0 - 1.013250 

5 2 

10 newtons/m * 1013.25mbar. This corresponds to 760mm Hg or approximately 
29.92 in Hg. In English units the pressure is 2116.22psf. 

Using this assumed temperature profile the pressure and density can be 
easily generated by integration. The analytical expressions are given by: 

p/p 0 - [ 1 + cth/T 0 ]‘ n h < llkm 

P/Pn - exp[ -(h - i )/h s 3 h > 11 km 
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where p 0 Is the assumed pressure at sea level (1013.25 mb) 
n - g/aR - 5.25588 

* RT/g - 6.3416km « 20,610ft (scale height) 

Pn“ 226.321 mb (Calculated pressure at 11 km) 
h lx ~ 11km - 36089ft (Tropopause altitude) 

One additional note is that in the model atmosphere the gravitational 

constant 'g' is taken to be constant in all equations and is defined as g * 
2 

9.80665 m/s . In fact, g varies with altitude (as well as with latitude), 
hence the resultant altitude in the model atmosphere Is referred to as 
'geopotential altitude' and is not the geometric (or physical) altitude. 
Since 'g' decreases with altitude, the geopotential altitude will be 
slightly less than the geometric altitude. The relation between geometric 
altitude and geopotential altitude is easily derived (See Ref. 1) but for 
the altitudes covered by commercial air traffic the difference between these 
two is considered insignificant. (At the tropopause the difference is about 
60ft.) In this report no distinction is made between geopotential and 
geometric altitudes. 

1.2 Atmospheric Variations 

Obviously there are few 'average 1 days where the temperature profile 
exactly matches the profile assumed In the model atmosphere. Also surface 
pressure Is highly variable from the standard value although percent 
deviations are smaller than those In temperature. Surface temperatures in 
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particular exhibit considerable variation from 15 deg C, but even at higher 
altitudes temperatures may deviate significantly from the model values. 
Figure 1, taken from Reference 1, depicts the historically observed 
statistical variations of atmospheric density and temperature. These 
fluctuations are large enough to have a significant effect on aerodynamic 
and thrust forces on an aircraft, hence consideration of an accurate ambient 
atmospheric is necessary for accurate trajectory prediction. 

The effect on aircraft performance is felt several ways. Since an 
aircraft altimeter is calibrated to translate pressure into an altitude 
reading based on the model atmosphere (see Section 1.4), the aircraft 
altimeter will not read the true geometric altitude. Nor will the 
calibrated airspeed measured in the cockpit, have the standard relationship 
to true airspeed due to an implicit temperature dependency. Engine 
performance will be different than expected since the air density and 
temperature are different from the standard values. Finally, wind errors 
will be introduced because of the difference between true altitude and the 
altimeter reading. What this translates to is that uncertainty in the 
temperature profile can be a major source of error in the prediction 
algorithms which determine expected landing times based on standard descent 
procedures . 

1.3 Incorporation of Actual Atmospheric Data into CTAS 

One of the primary elements of the CTAS (Center/Tracon Advisory System) 
software being developed by NASA is the trajectory prediction component 
which takes aircraft entering into the center advisory region and calculates 
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precise trajectories to touchdown. Key to the success of the CTAS system is 
its ability to project the arrival times at various waypoints so as to 
ensure smooth uninterrupted flow through the system. To ensure the accuracy 
of such a system two components are needed: 1) An accurate database of 
aircraft types, force coefficients and aerodynamic data, and 2) accurate 
atmospheric models with requisite spacial and temporal fidelity. One 
component of this atmospheric model is the accurate knowledge of pressure, 
temperature and density in the atmosphere. Another is the problem of 
determining accurate wind information which is discussed in the next 
section. 

From the previous section it is obvious that the atmospheric parameters 
are determined by the surface pressure and the variation of temperature with 
altitude. This latter requirement causes difficulty since often only 
surface temperatures are available. Standard correction procedures 
advocated for use by pilots use the standard lapse rate of the model 
atmosphere, coupled with surface temperature data, to compute density and 
true altitude variations from the pressure altitude, or altimeter reading. 
(See e.g. Ref. 2) This may also include either an adjustment of the 
tropopause height or an adjustment of the stratosphere temperature. Such 
corrections however may be significantly in error, since surface conditions 
are quite variable, and do not always reflect upper atmosphere changes. 

This is especially true for high altitude airports, such as Denver's 
Stapleton Airport which is above 5000ft. The model atmosphere temperature 
at this height is 5 deg C (41 deg F) while surface temperatures of 90 deg F 
are quite common in Summer. Projecting temperature variations at upper 
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altitudes based on this number alone introduces the possibility of 
significant errors. Even if surface pressure is 'standard', temperature 
deviations from standard values in the vicinity of the surface can introduce 
significant differences between the 'pressure altitude' (i.e. the altitude 
corresponding to a pressure level in the standard atmosphere) and the true 
altitude at much higher altitudes. Figure 2a shows a temperature profile 
which deviates from the standard in the altitude range from 0- 5000ft, then 
follows the model atmosphere temperature profile above 5000ft. Assuming the 
surface pressure is standard the deviation between the pressure altitude and 
the geometric altitude is shown on Figure 2(b). For this profile peak 
altitude deviations are about 1500 ft. 

If no high altitude temperature data is available, then some aprlori 
modeling procedure is clearly necessary. If however additional temperature 
data is available, as, for example, from descending or ascending aircraft, 
then these temperature variations can easily be included into an accurate 
atmospheric model . 

1,4 Implementation of an Atmosphere model for CTAS 

In this section we assume that a temperature profile and the surface 
pressure In the atmosphere are known. What we detail here is how to 
implement this into the CTAS algorithms. We assume that temperature is 
known at a series of altitudes (h^, T^), i— l,n. Then assuming a linear 

temperature variation between points, the pressure and temperature profiles 
are : 


9 



(T/T i ) - 1 + a.(h - h.)/T i 


p/ Pi - (T/T i ) ’® /a i R 
for h^ < h < h^ + ^and where 



,n-l 


At the last point h^, T^, or if * 0, (that is, if the atmosphere is 

locally isothermal) the pressure relation above is replaced by the 
exponential form 


p/ Pi - exp( -g(h - h i )/RT i ) 

Our procedure then is to sequential compute the temperature and pressure at 
the interface between successive layers, and to store the atmospheric 
exponent or the scale height for each layer. This array of numbers is 
computed by CURR_ATM (See Appendix 1) and is used by routines which 
specifically compute temperature and pressure as functions of altitude. 
Density calculation is not explicitly added but should be computed from the 
perfect gas law (p = p/RT) . 

If temperature data comes from aircraft aloft, then it should be 
recognized that the aircraft altitude will be a pressure altitude and not 
the geometric altitude required in the equations above. This can be 
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corrected by using the model atmosphere relations and the lower altitude 
temperature data to explicitly calculate the difference between geometric 
and pressure altitudes. 

The coding for the current atmosphere routines is done as a script file 
in Mat lab. The Mat lab code is 'C-like' however and a conversion to 'C' or 
some other higher order language is relatively straight-forward. See 
Appendix 1 for a source listing of the routines available, 

1.4.1 Note on Altimeter Settings 

An aircraft altimeter is calibrated to read out an altitude above mean 
sea level based on the model atmosphere. There is an additional setting 
parameter however which allows the altimeter to deviate from this value. 

For safety of flight, when operating below 18,000ft, pilots are required to 
adjust their altimeters to a fictitious pressure broadcast by the local 
station. This fictitious 'corrected' pressure is such that at the physical 
ground level the altimeter reads the actual height above sea level for the 
local point in question. For commercial air traffic this generally 
corresponds to the airport, from which the aircraft is departing from, or 
arriving at. This corrected barometric pressure is entered manually by an 
adjustment screw on the altimeter and displayed in the 'Kollsman window' of 
the altimeter. In the United States, the Kollsman window is calibrated to 
inches of mercury. Nominal sea level pressure is 29.92 in Hg and this 
setting is always used when flying above 18,000ft. In general the 
atmosphere Is not 'standard' and hence for an aircraft at the runway 
altitude, the altimeter will not read the posted runway altitude without 
some adjustment . The broadcast Kollsman window setting is simply a 
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fictitious pressure that rotates the altimeter dial so that a correctly 
calibrated altimeter reads correct runway height when the aircraft is on the 
runway. The physical interpretation of the Kollsman setting is shown in 
Figure 3, taken from Ref. 2. It is important to note that for a non- 
standard atmospheric pressure variation, the altimeter will read the correct 
geometric altitude only at the runway. At other altitudes the altitude must 
be corrected to obtain geometric altitude just as in the standard atmosphere 
case . 

2. Wind modeling 

Lack of accurate data on current atmospheric winds would be a serious 
impediment to the accurate predicting of aircraft arrival times at the 
runway and at various Intermediate waypoints. While this would not pose a 
safety problem for the CTAS system, the lack of accurate winds would result 
in increased controller advisories, and a loss of aircraft efficiency due to 
the extra maneuvering required to maintain spacing and arrival rates. 
Fortunately there are several sources of wind data available which make wind 
profile generation feasible. Wind data is available from: (1) 

Meteorological information available in the vicinity of most airports, (2) 
Enroute information from inertially equipped aircraft, and (3) Wind data 
from the NOAA Wind Profiler network which is available (or soon will be) for 
31 selected sites across the country. 

Access to, and usage of, the NOAA Profiler system data is discussed 
extensively in the next few sections of this report. Data processing 
requirements for the other data types (items (1) and (2) above) is similar; 
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only the methods of data collection are different. In the current study , 
data collection for these data types has not been considered and is a topic 
for further research. 

2.1 Wind data: NOAA Profiler system 


The wind profiler system currently being implemented by the National Oceanic 
and Atmospheric Administration (NOAA) was designed with the primary goal of 
weather prediction, particularly to help model and predict thunderstorm 
activity in the central section of the United States. There are currently 
31 sites planned. The locations of these sites are shown in Figure 4 and 
Table 2.1. For background on the Profiler system see Reference 3, 

The physical setup for a typical profiler site is depicted on Figure 5. 
At the Profiler site three doppler radar beams are used. One beam is 
vertical, the other two are each 15 degrees from the vertical and oriented 
orthogonal to each other, preferably in the North and East directions, but 
generally skewed due to potential interference with aircraft and satellite 
communications. (See Table 2.1). Back scattering of the radar beam from 
particulates and small scale turbulence in the atmosphere allow 
determination of the radial wind velocity along the three beams as a 
function of altitude. The altitude dependency is separated out by sending 
out radar pulses and processing reflected radar returns as a function of the 
time delay. These three non- orthogonal velocity components can be 
transformed to give magnitude and heading of the horizontal component of the 
wind, which are available as outputs of the Profiler station, as well as a 
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vertical component (which is not available to external users.) For details 
on the operation of the Profiler, see Ref. 3. 

The wave length of the radar determines the altitude resolution of the 
beam and its maximum altitude capability. Higher frequencies (hence shorter 
wave lengths) give better resolution but have a lower usable altitude range 
due to the lack of smaller scale turbulence at higher altitudes. The 
standard profiler site uses a frequency signal of 404 MHz (A « 70cm), which 
can give data, up to about 60,000ft. This system is operated in a high 
resolution mode and a low resolution mode by varying the pulse width and the 
averaging interval. In the low resolution mode, winds can only be obtained 
up to about 35,000ft. There is an additional Profiler site at Denver 
Stapleton airport that is not part of the Profiler Network. While 
functioning quite similarly to the Profiler sites, this site is an 
experimental station and operates at a frequency of 915 MHz. Hence it 
provides more accurate data over a smaller altitude range, and as discussed 
later, is more sensitive to environmental conditions.). 

At the time of this writing, only the station at Platteville, Colorado 
Is available for data, although several other stations are currently 
operating in a test mode. The complete complement of 31 sites should be 
available by the end of 1991. 

For use in wind prediction in the ATC problem, some knowledge of the 
data characteristics of this system are needed to aid in utilization of the 
data. Several factors can influence the quality of the received data. The 
primary adverse factor is precipitation. Since falling precipitation can 
produce a large return signal, any precipitation can seriously affect the 
quality of the data, particularly due to spatial inhomogeneities in droplet 
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speeds and distribution. This effect is especially true at higher radar 
frequencies, hence the Denver experimental facility frequently is most often 
affected by precipitation. 

This effect is magnified if local inhomogeneities affect the different 
beams in different ways. In fact, local inhomogeneities of any type can 
cause discrepancies since the velocity transformations to compute horizontal 
velocity assume spatial homogeneity. Due to the beam deflection angles this 
means that at high altitudes, beam locations are several miles apart. The 
velocity calculation effectively requires differencing to remove the 
vertical wind component. Any spatial variation in vertical winds between 
the different beams will be translated into erroneous horizontal velocity 
components . 

A third factor affecting the Profiler can be local changes In 
reflectivity of any type. For example the operation of the doppler radar Is 
dependent on backscattering of the emitted beam from dust and density 
fluctuations in the air. The reflection is particularly sensitive to 
fluctuations whose length scale Is half a wavelength. Particularly at high 
altitudes, inherent atmospheric stability and lack of turbulence can make 
the reflected signal extremely weak. In these cases backscattering from 
other sources in the side lobes of the radar can cause spurious reading. 

At the radar site a number of separate averages is performed to 
determine wind speed components at each altitude. As part of the Profiler 
output the number of averages and the returned signal power Is given as 
output data also. If insufficient averages are available, or the data fails 
a consistency check, the value '-999' is given for wind magnitude to 
indicate a bad data point. While the number of averages and the power value 


15 


should give some indication of data quality, our experience has been that 
this is a rather unreliable indicator. A better choice appears to be to 
view the data graphically , and to visually scan for bad data points. Data 
currently available at NASA Ames is for the Platteville and Denver sites. 
Typical data received from these stations is shown in Table 2.2. 

Some Profiler stations will be equipped to report on additional surface 
conditions. In particular, the temperature and pressure are of interest for 
the ATC problem. 

2.2 Wind Profile estimation 

Typical wind versus altitude data is shown in Figure 6 for the 
Platteville and Denver reporting stations. This data needs to be further 
smoothed, and outlier points removed for a consistent trajectory prediction 
algorithm. Several types of assumptions on wind variation with altitude 
could be used. For our purposes we have assumed that the data is to be fit 
to a series of contiguous straight line segments of velocity versus 
altitude. The number of segments is variable and is an input from the 
analyst at the time the data is to be processed. Fitting is done through a 
weighted least squares algorithm, which is identical to the prediction phase 
of a Kalman filter. Data is fit independently for the North and East 
components of the wind. Note that since the Platteville site is aligned in 
the North-East direction, this is equivalent to smoothing the raw data 
directly. For other sites with a non-North-East orientation, estimation 
applied to the natural orientation of the site Is probably preferable to 
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North-East. Note also that while the wind components each vary linearly 
with altitude, the magnitude and wind heading will not. 

For a set of fixed altitude break points for the various segments the 
data is fit to the segments by the weighted least square algorithm, (see 
e.g. Bryson[3]) Generally for Profiler data, all data weights for all wind 
values are assumed to be the same. Outliers are generally explicitly removed 
from the data although a much reduced data weighting would have virtually 
the same effect. A few typical fits to the data are shown in Figure 7. If 
desired, the breakpoint altitudes can also be included as parameters in the 
data fit. This can lead to improvement in the data residuals, at the 
expense of more lengthy and possibly ill-conditioned computations. A 
typical fit with variable breakpoints is shown in Figure 8. Note while this 
case is well-behaved, the fitting process here is non-linear, and hence an 
iterative improvement process is required to produce a final fit. If 
initial points are chosen reasonably, In our experience this fitting process 
usually proceeds with little difficulty. With a poor initial guess however, 
this fit can, and has, diverged in some of our simulations. Since 
divergence in an operational scenario could be quite serious, some 
safeguards to prevent this would be necessary In an operational setup. 

The details of the fitting procedure follow the standard least squares 
procedure. The Matlab code which has been developed to perform the fit for 
the Profiler data is given in Appendix 2. Note that in this code, we have 
added separate outlier logic which searches the fitted data for points which 
are significantly distant from the fit and automatically marks such data for 
deletion. In practice , this section worked quite well, and removed only 
those points which visually we felt to be outliers. Occasionally, if data 
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was particularly noisy and data residuals were high, then some points which 
visually appeared to be outliers were not detected by this algorithm. In 
such cases however, the entire Profiler wind data was often suspect, and 
consideration of other data types may be necessary. 

2.3 Other Data types 

Other wind data types should be available in an operational system. 

This will be necessary if the appropriate spatial resolution of the winds is 
to be determined, and will be necessary at sites not close to a wind 
profiler site. Such data is most likely data from inertially equipped 
aircraft, which can determine winds from the combination of air and inertial 
data. This data is currently available to the airlines- what is needed is a 
mechanism for obtaining this data for purposes of Air Traffic Control. 

Data processing of additional data types is considered to be identical 
to that of the Profiler data. The only decision is to estimate the relative 
accuracy of the data so that the appropriate data weighting relative to the 
other data types can be established. 

2.4 Spatial & temporal variation 

Wind profiles can exhibit considerable spatial and temporal variation. 
This variation must be accounted for and captured in a wind estimation 
strategy for CTAS . For example, a plot of some typical wind profiles 
measured at close time intervals at the Denver and Platteville profiler 
sites is shown in Figure 9. Similarly, overplots of the Denver and 
Platteville winds reported at the same times demonstrate the spatial 
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variability from the two sites (which are less than 50 mile apart) . See 
Figure 10. Because of this, the wind models proposed must be dynamic. That 
is, it is anticipated that once sufficient data becomes available, a fully 
time varying, three dimensional wind profile should be the goal of this 
system. Looking at aircraft arriving on major jetways, the initial goal is 
to formulate separate wind profiles on major jetways. Temporal variations 
are included by appropriate filtering of past and current wind data. In 
particular, a Kalman filter structure with added process noise to fade out 
old data is a straight-forward and reliable way to ensure optimum usage of 
wind data. Specific filter design characteristics such as required process 
noise will need to be determined once real data is available. 

3. Conclusion 

This report summarizes the author's initial efforts toward formulation 
of a wind and atmospheric model for Implementation into CTAS . A basic 
approach and set of algorithms has been developed for the wind and 
atmospheric estimation problems. Additional data types can easily be 
Included in the estimation algorithm. Attention to the logistics of data 
collection is a significant problem to be addressed in the future to achieve 
the level of accuracy sought in the CTAS system. 
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Table 2.1 Profiler sites 


Table 2.2 Raw data received from Denver/Platteville profiler 


SITE: DENVER — horizontal low mode WIND profile 
DATE: 90/07/23 

TIME: 20:00:00 (UTC) at start of acquisition period 

PROFILES: 12 

TIME DOMAIN AVGS : 124 

SPECTRA: 10 

PULSE WIDTH: 1.33 MICRO-SECS 

PULSE REP PERIOD: 50.00 MICRO-SEC'S 

MAX HOR VEL: 51.05 M/S 

FIRST HEIGHT: 289.78 M/AGL 

DELTA HEIGHT: 193 METERS 

HEIGHTS (GATES) : 32 




SPEED 

DIRECT 

HEIGHT 



POWER 

(EAST) 


GATE 

M/S 

DGRS 

KM/MS L 

#E 

#N 

DB 



1 

.9 

306 

1.90 

11 

9 

65.2 

C 

i = 

2 

2.8 

83 

2.09 

9 

7 

77.5 

C 


3 

4.1 

39 

2.29 

11 

6 

71.0 

C 


4 

.9 

43 

2.48 

11 

11 

73.1 

C 


5 

1.2 

326 

2.67 

10 

10 

74.8 

C 


6 

3.3 

307 

2.87 

12 

10 

75.0 

c 

__ 

7 

3.0 

299 

3.06 

12 

10 

75.1 

c 


8 

6.0 

302 

3.25 

12 

11 

67.3 

c 


9 

9.5 

281 

3.45 

12 

11 

59.8 

C 


10 

9.5 

290 

3.64 

12 

12 

57.0 

C 

— 

11 

9.9 

290 

3.83 

12 

12 

58.4 

c 


12 

11.3 

290 

4.03 

12 

12 

54.5 

c 


13 

11.9 

298 

4.22 

12 

12 

50.3 

C 


14 

12.4 

296 

4.41 

12 

11 

48.4 

c 

mmm 

15 

11.6 

300 

4.61 

12 

12 

55.3 

c 


16 

12.3 

304 

4.80 

12 

12 

56.9 

c 

I | 

17 

12.3 

309 

4.99 

12 

12 

52.6 

c 

j ? 

18 

12.8 

309 

5.18 

12 

12 

52.0 

c 


19 

13.4 

308 

5.38 

12 

12 

49.2 

c 


20 

13.7 

304 

5.57 

12 

12 

50.2 

c 


21 

14.1 

301 

5.76 

12 

12 

49.9 

c 


22 

14.6 

296 

5.96 

12 

11 

50.5 

c 


23 

14.1 

297 

6.15 

12 

10 

46.9 

c 


24 

14.7 

299 

6.34 

12 

10 

48.0 

c 


25 

15.7 

309 

6.54 

11 

11 

51.5 

c 


26 

17.1 

306 

6.73 

12 

12 

45.7 

C 


27 

18.1 

302 

6.92 

12 

12 

43.0 

C 


28 

18.3 

304 

7.12 

12 

12 

37.1 

c 

- - 

29 

18.5 

304 

7.31 

11 

12 

37.8 

c 

mm 

30 

18.5 

309 

7.50 

11 

12 

41.9 

c 


31 

18.4 

310 

7.70 

12 

12 

43.4 

c 


32 

18.1 

315 

7.89 

11 

12 

41.1 

C 


SITE: PLATTEVILLE — horizontal low mode WIND profile 
DATE: 90/07/23 

TIME: 20:00:00 (UTC) at start of acquisition period 

PROFILES: 12 

TIME DOMAIN AVGS: 350 

SPECTRA: 8 

PULSE WIDTH: 3.67 MICRO-SECS 

PULSE REP PERIOD: 238.00 MICRO-SECS 

MAX HOR VEL: 69.81 M/S 
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— 

Table 

2 . 2 Raw data 

received from Denver/Platteville profiler (cont) 

- ' 

FIRST 

HEIGHT 

i . 

1786.96 

M/AGL 



DELTA 

HEIGHT 

1 : 

289 

METERS 



HEIGHTS (GATES) : 

24 






SPEED DIRECT 

HEIGHT 



POWER (EAST) 


GATE 

M/S 

DGRS 

KM/MS L 

#E 

#N 

DB 


1 

2.8 

91 

3.31 

10 

12 

61.1 


2 

1.4 

323 

3.60 

11 

11 

61.5 


3 

4.7 

290 

3.89 

12 

11 

58.2 


4 

8.8 

297 

4,18 

12 

11 

62.6 


5 

11.1 

298 

4.47 

12 

11 

65.7 

— ‘ 

6 

12.2 

296 

4.76 

12 

12 

66.8 


7 

12.8 

295 

5.05 

12 

12 

66.1 


8 

13.1 

296 

5.34 

12 

12 

65.2 


9 

13.3 

298 

5.63 

12 

11 

66.4 

— 

10 

13.4 

303 

5.92 

12 

11 

68.5 


11 

13.4 

304 

6.21 

12 

11 

67.5 


12 

13.2 

306 

6.50 

12 

12 

62.4 


13 

14.3 

305 

6.79 

12 

12 

55.9 


14 

16.3 

302 

7.08 

12 

11 

50.2 


15 

15.4 

307 

7.37 

12 

12 

47.6 


16 

15.4 

313 

7.66 

12 

12 

50.1 


17 

16.2 

314 

7.95 

12 

10 

49.1 


18 

18.9 

309 

8.24 

11 

10 

47.5 


19 

19.3 

308 

8.53 

11 

8 

44.9 


20 

20.9 

306 

8.82 

11 

8 

40.3 

— 

21 

19.8 

300 

9.11 

12 

11 

38.9 


22 

20.1 

293 

9.40 

12 

12 

40.5 


23 

20.6 

290 

9.69 

12 

11 

39.5 

I 1 

24 

21.1 

292 

9.98 

12 

10 

33.7 


SITE: 

PLATTEVILLE 

— horizontal 

high mode WIND profile 

r : 

DATE: 

90/07/23 





mmm 

TIME: 

20:00: 

00 (UTC) at start of 

acquisition period 


PROFILES: 


12 




i i 

TIME : 

DOMAIN 

AVGS : 

124 




%} 

SPECTRA: 


16 





PULSE 

WIDTH: 


9.67 

MICRO-SECS 


PULSE 

REP PERIOD: 

672.00 

MICRO-SECS 


MAX HOR VEL: 


69.78 

M/S 



M 

FIRST 

HEIGHT 

: 

4201.77 

M/AGL 



DELTA 

HEIGHT 

z 

869 

METERS 



HEIGHTS (GATES) : 

16 




w 

SPEED DIRECT 

HEIGHT 



POWER (EAST) 


GATE 

M/S 

DGRS 

KM/MSL 

#E 

#N 

DB 


1 

13.3 

302 

5.73 

12 

11 

71.8 


2 

13.7 

304 

6.59 

12 

11 

67.8 


3 

15.0 

308 

7.46 

11 

11 

55.6 

= - 

4 

16.7 

311 

8.33 

11 

12 

51.0 

y 

5 

18.9 

298 

9.20 

12 

11 

44.6 


6 

20.7 

294 

10.07 

11 

10 

40.2 


7 

22 . 8 

292 

10.94 

10 

8 

35.6 


8 

23.2 

293 

11.81 

11 

10 

36.5 

— 

9 

23.7 

285 

12.68 

11 

9 

39.6 


10 

21.8 

285 

13.55 

11 

9 

39.8 


11 

19.5 

282 

14.42 

12 

10 

33.9 

*** 








— 










Table 2.2 Raw data received from Denver/Platteville prof Her (cont) 


12 

13.2 

303 

15.29 

11 

7 

31.0 

13 

9.3 

321 

16.16 

7 

11 

29.9 

14 

4.3 

261 

17.03 

7 

10 

26.6 

15 

6.6 

279 

17.90 

4 

5 

25.5 

16 

-999.0 

-999 

18.77 

5 

3 

27.5 






Figure 3 Explanation of the Kollsman window setting. 


WIND PROFILER SITES 



Figure 4 Profiler sites across the United States. 
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Figure 7 Typical straight line fits to Profiler wind data 
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Figure 9 Temporal variation of winds at Denver 


Altitude - kft Altitude - 


Comparison of Platteville(high) winds 



V EAST - kts 


V NORTH - kts 

















* II I 


APPENDIX A 


ZNASA Ames atmospheric routines (Matlab code) 
% contains: 


curr_atm.m 
inter .m 
pres_alt.m 
pressure . m 
rd_alt . m 
ttempf .m 
value . m 


Generates atmosphere from temperature data 
Interpolation utility 
Return pressure altitude given h 
Return pressure (uses curr_atm) given h 
Returns altimeter reading given pressure 
Returns temperature given h 
Interpolation utility 


% curr_atm.m gls 8/29/90 

ZPrograra to generate a matrix of current atmosphere parameters 
Z to determine pressure, density, temperature as functions of 
geopotential altitude. Output is matrix ATM defined by 


ATM— [ temp f ; 
alt; 
presO ; 
alph; 
nexp ] ; 


ZInput row vector of temperature (deg R) 

ZInput row vector of altitude (ft) 

ZOutput row vector of pressure (presO(l) input) 

ZOutput lapse rate between alt(i) -alt(i+l)/Hs at last pt 
ZPressure exponent or scale height if alph—0 


Used by function Pressure 


The current temperature, altitude data points are fixed here. 
Change to read from a file if desired 
ZN.B. HERE input tempf is deg F- -output is dee R 
tempf— [ 100 80 -69.7]; 
alt-[0 5000 36089] ; 
p0=1013 . 25 ; 


ZTHIS IS THE STANDARD ATMOSPHERE 

Ztempf=[59 41.169 -69.7]; Ztemp deg F for altitudes in h 
Zalt- 1000* [0 5 36.089]; Zaltitude in thousands of ft 
Zp0= 1013.25; % pressure at hO 

Z INTERNAL variables are all metric for easy reference 

cm2f - 1/.3048; Zmeter to ft EXACT 

gO - 9.80665; R - 1000*8.31432/28.9644; 

conv- 1 . 8/ cm2f ; Z deg K/m to deg R/ft 

gqR - conv*g0/R ; Z units now deg R/ft 

cN2p =1/(0 . 45359237*g0) ; ZNewtons to pounds 

tref- 1.8*273.15 -32; Z 0 deg F 

eps - l.e-7; Ztolerance for "0" lapse rate 

% end setup 


LL- length( tempf ) ; 

Z Form an array of lapse rates, exponents and pressures 
Z assume isothermal above last temperature eiven 
presO(l) = P 0; & 

for ii=l:LL-l 

alph(ii) - (tempf (ii+1) - tempf(ii) )/( alt(ii+l) - alt(ii) )• 
if abs( alph(ii)) > eps; 


A- 1 


nexp(ii) - -gqR/alph(ii) ; 

presO(ii+l) - presO(ii)*( (tref+tempf (ii+1) )/(tref+tempf (il) ) ) A nexp(ii) ; 
else 

nexp(ii) — (tempf(ii) +tref)/gqR; Xthis is scale height 
presO(ii+l) - presO(ii)*exp( - ( alt(ii+l) - alt (ii) )/nexp(ii) ); 
end; 
end; %for 
alph(LL)- 0; 

nexp(LL) - ( tempf (LL) +tref)/gqR; Xthis is scale height 
for ii-l:LL, tempf(ii)- tempf (ii)+tref;end; 

ATM- [ tempf ; al t ; presO ; alph ; nexp ] ; 

***************************************************************************** 

%function to return interpolation indices for one dimensional array 
% See function value to return interpolated value 
%function zzz — inter (h,harray) 

X returns zzz- [iref, a] such that 
X h - (l-a)*harray(iref) + a*harray (iref+1) ; 

% iref - 0 if h < harray(l) 

% iref - -1 if h > harray(max) 

% and a- h - harray(max) 

% gls 8/23/90 

function zzz - inter(h,harray) 

[nr ,nc]-size(harray) ; 
if min( [nr ,nc] ) —1 
Ll-max( [nr , nc] ) ; 

else disp('Not a 1-d array '); return 
end; 

if h>-harray (LI) 

iref-Ll; a- h-harray(Ll) ; 
elseif h<harray(l) , 

iref-0; a- h- harray(l) ; 

else 

stop-0; i— 2 ; 
while stop— 0; 
if hcharray(i) 
iref-i-1; 

a - (h - harray(iref ) )/(harray (iref+1) -harray(iref ) ); 
stop-1 ; 
else i-i+1; 
end; 
end; 

end; 

zzz-[ iref, a] ; 


***************************************************************************** 

% function to return pressure altitude for a given pressure 
X 

X altpres - pres_alt (press) ; OR pres_alt (press , str) str-'mb' or 'hg' 


A- 2 



X uses English units (ft, psf) unless mb - milibars 

J or hg - in hg specified 

% INTERNAL variables are all metric for easy reference 
X Use the NASA (ICAO) standard atmosphere 
X function altpres - pres_alt(press , str) ; 

function altpres - pres_alt(press , str) ; 

cm2f - 1/.3048; %meter to ft EXACT 

gO - 9.80665; R - 1000*8.31432/28.9644; 

gqR - gO/R; 

alphstd- -.0065; 

htropo - 11000; 

temps td- 288.15; %15deg C 

t_tropo- tempstd + alphstd*htropo ; 

hscale - R*t_tropo/g0 ; 

pstd- 1013.25; Xmillibars 

cN2p -l/(0.45359237*g0) ; XNewtons to pounds 

pstdE - 100*pstd*cN2p/(cm2f A 2) ; 

XpstdE— 2116.2; XEnglish units (psf) 

M2E- pstdE/pstd; 

p_tropo - pstd*(t_tropo/tempstd) A ( -gqR/alphstd) ; 

X end setup 
if nargin —2 

if str— 'mb', disp('mb') 

elseif str— 'hg', press - press*pstd/29 . 9213 ;disp( 'hg' ) 
end; 

else press - press/M2E; 
end; 

if press > p_tropo 

altpres — (tempstd/alphstd)*( (press/pstd) A (-alphstd/gqR) - 1); 
else altpres - htropo + hscale*log(p_tropo/press) ; 
end; 

altpres=cm2f*altpres ; Xremove this to return in meters 


% Given a temperature profile, generate pressure vs. altitude 
X Assume linear temperature profile between points 

X function press - pressure (alt , ATM) 

X alt in "ft"; pressure in "mb" 

X Matrix ATM generated by curr_atm.m contains current atmosphere model 

X Uses functions inter. m and value. m to get table values 
function press= pressure (alt , ATM) 

temp0= ATM(1 , : ) ; 
harray-ATM ( 2 , : ) ; 
pressO=ATM(3 , :) ; 
alph=ATM(4 , : ) ; 


A- 3 


nexp-ATM(5 , : ) ; 

zzz- inter (alt ,harray) ; 
iref-zzz (1) ; a-zzz (2) ; 
temph- value (tempO, zzz) ; 
if alph(iref) — 0 

press- pressO(iref )*( temph/tempO(iref ) ) A nexp(iref ) ; 
else press - pressO(iref)*exp( - (alt - harray(iref ) )/nexp(iref ) ) ; 
end ; 

***************************************************************************** 


% function to return altimeter reading for a given pressure 
X 

% altpres - td_alt(press , kollsman) ;0R rd alt(press , kollsman, str) str-'mb' or 
'hg' 

X press is the current pressure 

% English units (psf) unless mb - milibars 

* or hg - in hg specified 

X kollsman is the altimeter setting (kollsman window) 

X assumed in in Hg unless str-'mb' --> then in mb also 

X altpres is in (ft) 

X 

X Use function pres_alt.m which uses 
X the NASA (ICAO) standard atmosphere 
Xfunction altpres - rd_alt (press , kollsman, str) ; 

function altpres - rd_alt (press , kollsman, str) ; 

X function 
if nargln —3 

if str— 'mb' , 

hl-pres_alt (press , 'mb' ) ; 
h2-pres__alt(kollsman, 'mb' ) ; 
altpres- hi -h2; 
elseif str— 'hg' , 

hl*pres_alt(press , 'hg' ) ; 
h2-pres_alt(kollsman, 'hg' ) ; 
altpres- hi -h2; 
end; 

elseif nargin=l, kollsman=29 . 92 ; 
end; 

if nargin<3, 

hi- pres_alt(press) ; 
h2- pres_alt(kollsman, 'hg' ) ; 
altpres=hl-h2 ; 
end; 


***************************************************************************** 
% Given a temperature profile return temperature for a given h 


A-4 


% Assume linear temperature profile between points 
% 

X function temp-ttempf (h, ATM) 

X h is alt in "ft" 

X Matrix ATM generated by curr_atm.m contains current atmosphere model 

X Uses functions inter. tn and value. m to get table values 

function ttemp- ttempf (h,ATM) 

tempO- ATM ( 1 , : ) ; 

harray-ATM(2, :) ; 

zzz- inter (h .harray) ; 

iref-zzz(l); a-zzz(2); 

ttemp- value (tempO , zzz) ; 


****************************************************************-k-k-ki,ir******** 

X function to return function value using linear interpolation 
X on one dimensional array. 

X Input are interpolation indices 

X See function value to return interpolated value 

X 

X function h- value (harray, [iref, a]) 

X 

X [iref, a] returned from inter. m such that 
X h - (l-a)*harray(iref ) + a*harray(iref+l) ; 

X iref - 0 if h < harray(l) 

X returns h-harray(l) 

X iref — max if h > harray(max) 

X returns h-harray(max) 

X gls 8/23/90 

function h — value (harray, zz) 
iref- zz(l) ; a- zz(2); 

[nr ,nc]— size (harray) ; 
if min([nr,nc]) —1 
LI -max ( [ nr , nc ] ) ; 

else disp('Not a 1-d array' ) ;re turn 
end; 

if iref<=0 

h-harray(l) ; 
elseif iref—Ll 
h= harray(iref) ; 
else 

h - (l-a)*harray(iref) + a*harray(iref+l) ; 
end; 


***************************************************************************** 


A-5 


APPENDIX B 


XNASA Ames routines for wind profile fit (Matlab code) 

X contains : 

wnd7 .m : Read Profiler data and plot 
fit7.m : Least squares fit to data 

both these current versions are hardwired 
to eliminate options available in early code 

************************************************************** + ***^^^^ 

X wnd7.m 8/22/90 Master in /matlab/gls ALL Plat mod 
X 

Xtest wind profile 

Xload wind profile from data file 

disp( 'Data should be loaded in standard format') 

PLAT- 'PLAT'; DENV-'DENV'; 

El-exist( ' txtla' ) ; wherl-' '; 

E2— exist( ' txt2a' ) ; wher2— ' '; 

E3-exist( ' txt3a' ) ; wher3-' '; 
nam-[ 'V EAST ' 

'V NORTH' 

'V MAG ' 

'V E/N ' ] ; 

disp(' Plot data for') 

if El, wherl - txtla(7:10); disp([' 1: ' txtla] );end; 

if E2 , wher2 - txt2a(7:10); disp([' 2: ' txt2a]);end; 

if E3 , wher3 - txt3a(7:10); disp([' 3: ' txt3a]);end; 

disp([' 4: ' 'Platt.- ALL'j);end; 
disp( ' ') 

numin- input (' -->'); 

if numin—4, %this finds the Plattville data- -IF IT EXISTS- - 
platl=0 ;plat2«0 ; 
allplat-1 ; 

if wherl—PLAT, plat 1-1 ; numin- 1 ; 
end; 

if wher2— PLAT, 

if platl— 1, plat2-2; 
else platl~2 ; numin=2 ; 
end; 

end; 

if wher3 — PLAT, 

if platl>0 , plat2=3; 
else platl=3 ;numin=3 ; 
end; 

end ; 

else allplat=0; 
end; %numin-=4 

data= eval ( ['data' num2str (numin) ] ) ; 

txta- eval ( [ 9 txt 9 num2str (numin) 'a']); 

txta=[ txta(7 : 15) ' ' txtb(l:15) # ' txtc(6:14)]; 
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if allplat, 

[lowlen, temp] - size(data) ;%Save this for future use in fit 
data- [data;eval( ['data' num2str (plat2) ] ) ] ; 
txta- [ txta ( 1 : 1 5 ) '-ALL' txta(16 : 35) ] ; 
numin -4; 

[hilen, temp] -size (data) ; 

end; 

disp(' ') 

% EAST & NORTH winds only 
iplt-0; 

vmps- data( : , 2) ; 
degs- data(: , 3) ; 
altkm-data( : ,4) ; 


% 

cd2r« pi/180. ; 
cm2f~ 39.2/12; 
cfps2kt- 1/1.69; 

% make invalid data to nan 
LL - length (vmps) 
for ii«l:LL 

if vmps(ii) < 0, 

degs(ii)-nan; vmps(ii)- nan; 
end; 

end; %ii loop 

vxfps= -cm2f*vmps . *sin(cd2r*degs) ; 
vyfps- -cm2f^vmps . *cos(cd2r*degs) ; 
altkft« cm2f*altkm; 
vplot* cfps2kt*vxfps ; 
vplot2- cfps2kt*vyfps ; 

%make sure hold is off and screen is clear 

hold off 

clg 

subplot (121) ; % Do multiple plots 

if allplat, 

plot (vplot , altkft , '+' ) ;hold 
plot(vplot(l: lowlen) , altkft (1 ; lowlen) , ' : # ) ; 
plot (vplot (lowlen+1 : hilen) , altkft (lowlen+1 :hilen) / : ' ) ; 
else plot (vplot, altkft, ' : ' , vplot , altkft, '+' ) ; 
end; 

title(txta) ;grid; 

xlabel([nam(l, :) ' - kts']); ylabel ( 'Altitude - kft'); 

If taxis— axis ; axis ; Xcapture axes then release 
hold off 
if allplat, 

subplot(122) 

plot(vplot2 , altkft , '+' ) ;hold; 
plot(vplot2(l: lowlen) , altkft (1 : lowlen) ,':'); 
plo t(vplot2( lowlen+1 : hilen) .altkft (lowlen+1 rhilen) ,':'); 
else plot(vplot2 , altkft , ' : ' , vplot2 , altkft , '+' ) ; 
end; 

xlabel ( [nam(2 , : ) ' - kts']); 
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grid; 

rgtaxis-axis ; axis ; Xcapture axes then release 
X[xc,yc]=ginput(l, 'sc') ; 

%text(xc ,yc , txtb , 'sc') 

%text (xc , yc+ . 05 , txtc , ' sc ' ) 
hold on 

disp('do auto fit') 
fit7 


Xfit7.m 

% liLZ l ™i d t t& r t0 b ? S f S , trai ^ ht line se § n,ents using least square fit 
Break points for altitudes are variable in this version 

and number of segments are fixed. s 

See fitl.m for fixed break points. 


8/19/90 mod for fit7.m in iter— =1 only remove outliers automatically 
i /■ /i « y,L fitted data point is very far out of bounds 

f / i'o^oa 8 ^j 2 _ mod . to fit4 7/18/90 for fixed or variable 
6/19/90 add for fit on subplots -->fit5.m 

6/20/90 add weighted Is + add other vel . data 
keyboard command disabled due to problems on SUN's 


sigls-1; %Covariance of basic wind data 
sig2s-=. 01 ; XCovariance of "added” data 
epsp- 1 . e - 3 ; parmchg-1 ; 
hsegl=[ 25.0]; Xaltitudes in kft 
hMAX- 1000; 


dispJ'This program will fit straight lines to standard data' 
disp( Assume equal noise covariance of siels-'^ 
disp(sigls) 6 ' 


disp( ' It will also pick up additional data 
disp( ' VEAST WORTH HNEW if EXISTNEW is a 
disp('This data has covariance sie2s“') 
disp(sig2s) 

disp( 'Assume for this case data is plotted 
itermax-=10 Xfixed for now 
disp( ' itermax=l for fixed breakpoints'); 
disp ( ' ' ) 

itermax-input( 9 itermax- ') 


from arrays 9 ) 
defined variable' 


in E/N form') 


disp( 'Determine least squares fit') 

disp( ' to straight line wind profile segments') 

hsegl=input(' Input initial breakpoint vector- ')• 

nsegs- length(hsegl)+l ; 

nparm= nsegs+1 ; 

NEWDATA - exist( 'EXISTNEW'); XNEWDATA-1 if data 
if NEWDATA, lnewdata- length(VEAST) ;end; 

Xnumber of parameters (per curve) in ’ fit ’ (exclude 
if iplt—=0 | iplt— 4, nparm2=nparm; 

else nparm2-0;end; 
nparmt= nparm+nparm2 ; 


breakpoints ! 


XGet rid of nan's in data (+ outliers marked as nan) 


) 

) 


) 
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LL- vp lot— nan ; 
vf it-vplot (LL) ; 
h- altkft(LL); 

% Make sure its a column vector 

[nr ,nc]-size(vfit) ; if nr— 1 vfit- vfit'jend; 

[nr ,nc]-size(h) ; if nr— 1 h * h';end; 
if NEWDATA, vf it-[vf it ;VEAST] ; h-[h;HNEW]; end; 

%if subplots, alt data string is same on both 
if nparm2 > 0, vfit2- vplot2(LL) ;h2 - h; 
if NEWDATA, vfit2-[vf it2 ; VNORTH] ; end; 
vfitt-[vfit ; vfit2] ;ht«[h;h2] ; 
else vfitt-vfit; ht-h; end; 

% 

onevec-ones(vfit) ; 
onevect- ones(vfitt); 

% now work with only reduced data 
lvel- length(vfit) ; 
lold- lvel- lnewdata; 

% 

iter - 1; 

X 

X Enter main loopsize( 

% 

while parmchg > epsp %loop over estimator until convergence is met 
hseg- [ 0 hsegl hMAX] ; ^segment boundaries 
A- zeros (lvel ,nparm) ; 

A( : , l)-ones(vfit) ; 

L0= h > hseg(l)*onevec ; 
for jcol=2 :nsegs+l 

LI- h > hseg(jcol)*onevec ; 

L01 - L0 - LI; 

A(:,jcol) - hseg(jcol)*Ll - hseg(jcol-l)*L0 + h.*L01 ; 

L0 - LI; 
end; %for jcol 
if nparm2 > 0, 
zA- 0*A ; 

A= [A zA; zA A] ;ht-[h;h] ; 
else ht-h; end; 

X Now add columns for breakponts 
% but only in the second and higher iteration 
if iter —1, oldparm=zeros(nparmt,l) ; 
else %for iterations above 1, estimate break points 
A12- zeros(lvel ,nsegs-l) ; A22- A12; 
for j col-2 :nsegs 

LI- h > hseg(jcol)*onevec ; 

A12( : , jcol-1) - ( parm(jcol)- parm(jcol+l) )*L1; 
end; %for jcol 
if nparm2 > 0 

for jcol=2:nsegs 

LI- h > hseg(jcol)*onevec ; 

A22( jcol-1) - (parm(nparm +jcol)- parm(nparm+j col+1) )*L1; 
end; %for jcol 
A12- [A12 ; A22 ] ; 
end; %nparm2 
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A- [A A12 ] ; 

if iter— 2, oldparm=[parm; hseg(2 :nsegs) ' ] ; 

else oldparm - parm; 
end; 

end; Xif iter 

Vfit— A( : , 1 :nparmt)*oldparm(l :nparmt) ; ^Estimated measurements 

RIA-A; 

if NEWDATA, 

for irow-l:lold, 

RIA(irow,:)- RIA(irow, : )/sigls; 

RIA(lvel+irow, : )- RIA(lvel+irow, :)/sigls; 
end; 

for irow-lold+1 : lvel ; 

RIA( irow, : )- RIA(irow, : )/sig2s ; 

RIA(lvel+irow, :)- RIA(lvel+irow, :)/sig2s; 
end; 

end; Xif NEWDATA 
Api- inv(A'*RIA)*RIA' ; 
dparm-Api*( vfitt - Vfit); 
parmchg- sqrt(dparm'*dparm) ; 
parm- dparm+oldparm; 

XGet new residuals for test 

Vfit- A( : , 1 :nparmt)*parm(l :nparmt) ; XEstimated measurements 

delv - ( vfitt - Vfit) ; 

rmsdv - sqrt( delv'*delv/length(delv) ) ; 

maxdelv- max( abs(delv) ); 

disp(['RMS fit deviation - ' num2str( rmsdv) ' kts']) 
disp(['MAX fit deviation - ' num2str (maxdelv) ' kts']) 

X Check velocity deviations here on first iteration 
X auto delete outliers 
if iter — 1 

rmst— max([rmsdv 2.5]); XUse 4*rms or 10 as outlier test 
Lol- abs(delv) > 4*rmst*onevect ; 

lLol- Lol'*Lol; Xthis dot product gives the number of outliers 
if lLol — 0, disp('N0 OUTLIERS FOUND') 
else 

disp( [ num2str(lLol) ' OUTLIERS LOCATED AT alt-']) 
disp(ht(Lol) ) 

Xautodelete outlier- - tricky logic as we will have to delete 
Xdata from both sets if E/N data being used 
if nparm2 > 0 

Loll- Lol(l : lvel) ; Lol2- Lol(lvel+l : 2*lvel) ; 

Loll = Loll + (Lol2 > Loll); 

Lol2 - Lol2 + (Loll > Lol2) ; 

Loll - onevec - Loll; 

Lol2 = onevec - Lol2 ; XHere Loll & Lol2 should be the same 
Lol = [Loll;Lol2] ; 

XNow get rid of outlier 
vfitt- vfitt(Lol) ; h-h(Loll); 
vfit-vfit(Loll) ; vfit2- vfit2(Lol2); 

else 

Lol - onevec - Lol ; 

vfitt- vfitt(Lol) ; 

vfit- vfit(Lol) ; h= h(Lol) ; 
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end; %nparm2 

iter - 0; Xfinal step- set iter back so iter-1 done again 
onevect - ones(vfitt); 
onevec-ones (vf it) ; 

lvel— length(vfit) ;X****N0TE ASSUME FOR NOW DELETED DATA WAS IN 

OLDDATA SET 

end; Xif lLol 
end; X if iter— 1 

% if iter > 0 disp ( 9 Estimated parameters:') 

% disp(parm) 

% end; 

if iter >1, hsegl- parm(nparmt+l :nparmt+nsegs-l) ' ; 
end; 

iter- iter+1; 
oldparm- parm; 

status- [ iter, parmchg, hsegl ] ; %display output 
disp(status) 

if iter >itermax 
disp( ' - -quit' ) ; 

parmchg- -parmchg; end; X if itermax 
end; %while loop 
% 

% converged or quit. Now plot fit 
X 

htt- max( h ) ; 

hplot«[0 hsegl htt]; %plot breakpoints and last point 
vfit-[parm(l) ] ; 
for ii-2:nparm 

vfit(ii)- vfit(ii-l) + parm(ii)*(hplot(ii) - hplot(ii-l) ); 
end; 
hold on 
if nparm2 > 0, 

subplot (121) ; axis (If taxis) ;plot(vfit ,hplot , 'b' ) 
vf it2-[parm(nparm+l) ] ; 
for ii-2:nparm 

vfit2(ii)- vfit2(ii-l) + parm(nparm+ii)*(hplot(ii) - hplot(ii-l) ); 
end; 

subplot(122) ; axis (rgtaxis) ; 
plot(vfit2 ,hplot , 'b' ) 
else plot (vf it ,hplot , 'b 9 ) 
end; 

disp('End of fit7.m') 
end % fit program 

************************************************************************** 
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